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Abstract 

We show that it can be considered some of Bach pitches series as a stochastic process with 
scaling behavior. Using multifractal deterend fluctuation analysis (MF-DFA) method, frequency 
series of Bach pitches have been analyzed. In this view we find same second moment exponents 
(after double profiling) in ranges (1.7 — 1.8) in his works. Comparing MF-DFA results of original 
series to those for shuffled and surrogate series we can distinguish multifractality due to long-range 
correlations and a broad probability density function. Finally we determine the scaling exponents 
and singularity spectrum. We conclude fat tail has more effect in its multifractality nature than 
long-range correlations. 

PACS numbers: 02.50.Fz, 05.45.Tp 
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Pythagoras knew it, but Bach demonstrated it: without mathematics there is no music 



INTRODUCTION 



Many people think that, mathematics and music have some vague sort of affinity, but 
most often supposed relationship between two fields turns out to be in details that are 
not central to either. The mathematical proportions in musical works are hidden from 
listener from old days. Thus it is forced to make use of interpretative techniques in order 
to search for them, which is problematic from a methodological point of view. Some music 
historians have very little time for numerology. In addition to mathematics being seen as 
numerical symbolism, music is closely linked to absolute physical entities, such as frequency 
and relation between intervals (an interval is space between two notes). Already in antiquity 
this was seen as natural or cosmic premise on which music relied. It is illustrated fact that 
not just musical notation, but also relationship between music and time has something to 
do with mathematics and with one of the most significant transformations in music history. 
Complexity of music can be especially attracting scientific interest. Among great variety 
of complex and disordered systems most of music parameters such as frequency and pitch 
(Pitch is the sound frequency of any given note.) Amplitude or Dynamics 

(Dynamics are the changes in volume during a musical piece.) [7|, [8|, Intervals (Intervals are 
the distances between notes in the musical scale.), Rhythm (Rhythm is the structure of the 
placement of notes in musical time.) can be consider as stochastic processes. Also, some 
authors try to clustering the music j^], \\\ . 

In all technicalities, music can be composed of notes. A note is a sign used in music 
to represent the relative duration and pitch of sound. In traditional music theory pitch 
classes are represented by the first seven letters of the Latin alphabet (A, B, C, D, E, F, 
and G) or [Do - Re - Mi - Fa - Sol - La - Si), in Italian notation. Each note is assigned 
a specific vertical position on a staff. Since the physical causes of music are vibrations of 
mechanical systems, their frequencies are often measured in hertz (Hz). These frequencies 
are mathematically related to each other, and are defined around the central note. The 
current "standard pitch" for this note is 440 Hz. Any note is exactly an integer number 
of half-steps away from central note. Let this distance be denoted n. Then, the desired 
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FIG. 1: Typical sheet music, Invention no.l by Bach. 



frequency is given by 



v 



440 x 2 n/12 Hz. 



(1) 



In this paper we would like to characterize complex behavior of frequency of note signal of 
Bach's Inventions and Sinfonias through computation of signal parameters scaling exponents, 
which quantifies correlation exponents and multifractality of signal. Inventions and Sinfonias 
are collection of short pieces which Bach wrote for musical education of Bach young pupils. 
These are among finest examples of artistic gems ever written for this purpose, and probably 
because of this, they became very popular among Bach's pupils and others ever since they 
were written. Inventions and Sinfonias contain two and three music voices respectively. The 
number of the data in frequency series is dependent to the piece and will obtain from one 
voice which contains 1000 notes in average. Because of non- stationary nature of frequency of 
music series, and due to fmiteness of available data sample, we should apply a methods which 
are insensitive to non-stationarities, like trends. In order to separate trends from correlations 
we need to eliminate trends in our frequency data. Several methods are used effectively for 



this purpose: detrended fluctuation analysis (DFA) rescaled range analysis (R/S) [3] 
and wavelet techniques (WT) [131 ] . 

We use MF-DFA method for analysis and eliminating trend from data set. This method 
are the modified version of DFA method to detect multifractal properties of time series. 
DFA method introduced by Peng et al. 11] has became a widely used technique for the 
determination of (mono-) fractal scaling pro perties and the detection of long-range correla- 



tions in noisy, non- stationary time series 



diverse fields such as DNA sequences 
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171 ] . It has success: 



181 ]. heart rate dynamics 19 



ull y b een applied to 
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211 ] . neuron spik- 
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221 ]. human gait [23|, wind speed, [241 ] . long-time weather records [25(, cloud structure 



26j | . geology |27|, ethnology [28(, economical time series 



291 ] . and solid state physics 



30]. 



One reason to employ DFA method is to avoid spurious detection of correlations that are 
artefacts of non-stationarity in time series. 

The focus of present paper is on intriguing statistical properties and multifractal nature 
of frequency series. In Particular, Figures 1 and 2 show the score of Invention no.l and the 
frequency fluctuation of Invention no. 1 and Sinfoina no. 1, respectively. In general, two 
different types of multifractality in frequency series can be distinguished: (i) Multifractality 
due to a fatness of probability density function (PDF) of time series. In this case mul- 
tifractality cannot be removed by shuffling the series, (ii) Multifractality due to different 
long-range correlations in small and large scale fluctuations. In this case data may have a 
PDF with finite moments, e.g. a Gaussian distribution. Thus corresponding shuffled time 
series will exhibit mono-fractal scaling, since all long-range correlations are destroyed by 
shuffling procedure. If both kinds of multifractality are present, shuffled series will show 
weaker multifractality than original series. 

The paper is organized as follows: In section II we describe MF-DFA methods in detail 
and show that, scaling exponents determined by MF-DFA method are identical to those 
obtained by standard multifractal formalism based on partition functions. In section III, 
we analysis of frequency series of Bach's Inventions also examine source of multifractality 
in frequency data by comparison MF-DFA results for remaining data set to those obtained 
by MF-DFA for shuffled and surrogate series, section IV closes with a conclusion. 



II. MULTIFRACTAL DETRENDED FLUCTUATION ANALYSIS 



The simplest type of multifractal analysis is based upon standard partition function multi- 
fractal formalism, which has been developed for multifractal characterization of normalized, 
stationary measurements 
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34] . Unfortunately, this standard formalism does not 



give us correct results for non- stationary time series that are affected by trends or those which 
cannot be normalized. In the early 1990s an improved multifractal formalism wavelet trans- 
form modulus maxima (WTMM) method [lj, has been developed. This method is based 
on wavelet analysis and involves tracing the maxima lines in continuous wavelet transform 
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FIG. 2: Up to down: typical frequency series of Invention and Sinfonias no.l by Bach. 



over all scales. Other method like, multifractal detrended fluctuation analysis (MF-DFA), is 
based on identification of scaling of gth-order moment depending on signal length, and it is 
generalization of standard DFA method in which q = 2. In contrast to WTMM method MF- 
DFA does not require modulus maxima procedure, and hence does not require more effort in 
programming and computing time than conventional DFA. On the other hand, one should 
find correct scaling behavior of fluctuations, from experimental data which often affected by 
non- stationary like trends. This have to be well distinguished from intrinsic fluctuations of 
system. In addition often in collected data we do not know reasons, or even worse scales, for 
underlying trends, and also, usually available record data is small. So for reliable detection 
of correlations, it is essential to distinguish trends for intrinsic fluctuations from collected 



data. Hurst rescaled-range analysis 



12| and other non-detrending methods work well when 



records are long and do not involve trends, otherwise it might give wrong results. DFA is a 
well established method for determining scaling behavior of noisy data, where data presence 



of trends and we don't know their origin and shape 11 



20 



5 



A. Description of MF-DFA method 



Modified multifractal DFA (MF-DFA) procedure consists of five steps. The first three 

nnnnn 

steps are essentially identical to conventional DFA procedure (see e.g. [11|, UM, UM, Il7|). 
Suppose that Xk is a series of length N, and it is of compact support, i.e. Xk = for an 
insignificant fraction of the values only. 

• Step 1: Determine the "profile" 

i 

y(i)=£ [a* -(*)], i = l,...,N. (2) 

fc=i 

Subtraction of the mean (x) from Xk is not compulsory, since it would be eliminated by later 
detrending in third step. 

• Step 2: Divide profile Y(i) into N s = int(A/s) non overlapping segments of equal lengths 
s. Since length N of series is often not a multiple of considered time scale s, a short part at 
the end of profile may remain. In order not to disregard this part of series, same procedure 
should be repeated starting from the opposite end. Thereby, 2N S segments are obtained 
altogether. 

• Step 3: Calculate local trend for each of 2N S segments by a least-square fit of series. Then 
determine the variance 

F 2 (s, {y[iy - 1)« + i] - Vu®} 2 . (3) 

s i=i 

for each segment v — 1 , . . . , N s and 

F 2 (s, = {YIN -(u-N.)8 + {[- yS)Y , (4) 

s i=X 

for v — N s + 1, . . . , 2N S . Where, y v {i) is fitted polynomial in segment v. Linear, quadratic, 
cubic, or higher order polynomials can be used in the fitting procedure (conventionally 



called DFA1, DFA2, DFA3, . . ., DFAm) ll|, |21|. Since detrending of time series is done by 



subtraction of fited polynomial from profile, different order DFA differ in their capability of 
eliminating trends in series. In (MF-)DFAm trend of order m in profile (and equivalently, 
order m — 1 in original series) are eliminated. Thus a comparison of results for different 

nn 

orders of DFA allows one to estimate type of the polynomial trend in time series |15l . Ilq . 
• Step 4'- Average over all segments to obtain g-th order fluctuation function, defined by: 



F q (s) 



Irj:r-( S ,^' 2 \ , (5) 



2N, 



where, in general, variable q can take any real value except zero, q = 2, standard DFA 
procedure is retrieved. Generally we are interested to know how generalized q dependent 
fluctuation functions F q (s) depend on time scale s for different values of q. Hence, we must 
repeat steps 2, 3 and 4 for several time scales s. It is apparent that F q (s) will increase with 
increasing s. Of course, F q (s) depends on DFA order m. By construction, F q (s) is only 
defined for s > m + 2. 

• Step 5: Determine scaling behavior of fluctuation functions by analyzing log-log plots of 
F q (s) versus s for each value of q. If series Xj are long-range power law correlated, then F q (s) 
, for large values of s, increases as a power-law i.e., 

F q (s) ~ s fc fo>. (6) 

In general, exponent h(q) may depend on q. For stationary time series such as, fractional 
Gaussian noise (fGn), Y(i) in eq. El will have a fractional Brownian motion (fBm) signal, 
so, < h(q = 2) < 1.0. The exponent h(2) is identical to well known Hurst exponent H 



31 |. Also for non- stationary signal, such as fBm noise, Y(i) in eq. [2], will be a sum 



of fBm signal, so corresponding scaling exponent of F q (s) is identified by h(q = 2) > 1.0 
[ill . [3^]. For monofractal time series, h(q) is independent of q, since scaling behavior of 
variance F 2 (s, v) is identical for all segments u, and averaging procedure in eq.© will just 
give this identical scaling behavior for all values of q. If we consider positive values of q, the 
segments v with large variance F 2 (s, v) (i.e. large deviation from the corresponding fit) will 
dominate average F q (s). Thus, for positive values of q, h(q) describes scaling behavior of 
segments with large fluctuations. For negative values of q, segments v with small variance 
F 2 (s, v) will dominate average F q (s). Hence, for negative values of q, h(q) describes scaling 
behavior of segments with small fluctuations. 



B. Relation to standard multifractal analysis 

For a stationary, normalized series multifractal scaling exponent h(q) defined in eq.([6]) 
are directly related to scaling exponent r(q) defined by standard partition function based 
on multifractal formalism as shown below. Suppose that series of length iV is a sta- 
tionary, normalized sequence, then detrending procedure in step 3 of MF-DFA method is 
not required, since no trend has to be eliminated. Thus, DFA can be replaced by standard 
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fluctuation analysis (FA), which is identical to DFA except definition of variance, which is 
simplified for each segment v — 1, . . . , N s . In step 3 eq.Q now becomes : 

Fl A {s,v) = [Y{vs)-Y{{v-l)s)f. (7) 

Inserting this simplified definition into eq.© and using eq.Q, we obtain 



2N S 

2iV, 



1/9 



±<£\Y(v8)-Y({u-l) a )\*\ (8) 

For simplicity we can assume that length N of series is an integer multiple of scale s, 
obtaining N s = N/ s and therefore 

N/s 

\ Y ("s) - Y((u - l)s)\ q ~ s^- 1 . (9) 



This corresponds to the multifractal formalism used e.g. in 



32 



34[. In fact, a hierarchy of 



exponents H q similar to our h(q) has been introduced based on eq.(|H]) in [32j. In order to 



relate this to standard textbook box counting formalism [31|, [33], we employ definition of 
profile in eq.Q. It is evident that term Y(ys) — Y((u — l)s) in eq.([9]) is identical to sum of 
numbers Xk within each segment v of size s. This sum is known as box probability p s {v) in 
standard multifractal formalism for normalized series Xk, 

us 

Vs {v)= x k = Y(us)-Y((y-l)s). (10) 

fe=(iy-l) S +l 

The scaling exponent r(q) is usually defined via partition function Z q (s), 

N/s 

where q is a real as in MF-DFA method, discussed above. Using eq.l fTOl) we see that eq.l fTTl) 
is identical to eq.([H]), and one can obtain analytical relation between two sets of multifractal 
scaling exponents, 

r{q) = qh(q) - 1. (12) 

Thus, we see that h(q) defined in eq.([6]) for MF-DFA is directly related to classical mul- 
tifractal scaling exponents r(q). Note that, h(q) is different from generalized multifractal 
dimension 

q — 1 ?—l 
8 




FIG. 3: The log-log plot F(s) versus s for q = 2.0 after double profiling for Invention no. 1. 

that are used instead of r(q) in some papers. In this case, while h(q) is independent of q for a 
monofractal time series, D(q) depends on q. Another way to characterize a multifractal series 
is looking to singularity spectrum f(a), which is related to r(q) via a Legendre transform 

HQ, 

a = r'(q) and f(a) = qa — r(q). (14) 

Here, a is the singularity strength or Holder exponent, while f{a) denotes dimension of 
subset of series that is characterized by a. Using eq.(|T2l). we can directly relate a and f(a) 
to %), 

a = h(q) + qh' (q) and /(a) = q[a — h(q)] + 1. (15) 

Holder exponent denotes monofractality, while in multifractal case, different parts of struc- 
ture are characterized by different values of a, leading to existence of spectrum f(a). 



III. ANALYSIS OF MUSIC FREQUENCY SERIES 

As mentioned in section II, a spurious of correlations may be detected if time series is non- 
stationarity, so direct calculation of correlation behavior, spectral density exponent, fractal 
dimensions etc., don't give us a reliable results. It can be checked out that, frequency series 
is non-stationary. One can verified non-stationarity property experimentally by measuring 
stability of average and variance in a moving window for example using scale s. According 
to MF-DFAl method, generalized Hurst exponents h(q) in eq.([6]) can be found by analyzing 

9 
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FIG. 4: The q dependence of the exponents h*(q), T*(q) and singularity spectrum /*(«*), after 
double profiling, are shown in the upper to lower panels respectively for Invention no. 1 . 

log- log plots of F g (s) versus s for each q. Hurst exponent is between < H < 1. However, 
the MF-DFA method can only determine positive generalized Hurst exponents, in order to 
refine the analysis near the fGn/fBm boundary or strongly anti-correlated signals when is 
close to zero. The simplest way to analyze such data is to integrate the time series before the 
MF-DFA procedure. Hence, we replace the single summation in Eq. (1), which is describing 
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FIG. 5: The second moment exponents, h*(q = 2), after double profiling, of various Inventions 
and Sinfonias in their number ordering. 

the profile from the original data, by a double summation by using Classification by the 
signal summation conversion method (SSC). After using SSC method, fGn switch to fBm 
and fBm switch to sum-fBm. In this case the relation between the new exponent, h*(q = 2), 



and h(q = 2) is h(q = 2) = h*(q = 2) - 1 



38 



(recently Movahed et al. have proven 



the relation between derived exponent from double profile of series in DFA method and 
h(q = 2) exponent in Appendix of 39|). We find h*(q = 2) = 1.74 ± 0.03 of Invention no. 
1 by using SSC method. Note that, the exponents of the new series, after double profiling 
(SSC), are different from before SSC. Therefore we name these new exponents as h*(q), 
r*(q) and /*(«*). 

Results using MF-DFAl method for frequency signal are shown in Figure HJ which shows 
that frequency series is a multifractal process as indicated by strong q dependence of the 
exponents h*(q = 2) and T*(q) 40|. q dependence of the multifractal scaling exponent T*(q) 
has linear behaviors for q > and q < and slopes of T*(q) are 1.53 ± 0.01 and 2.61 ± 0.01, 
respectively. The values which derived for quantities of MF-DFAl method for Invention no. 
1 are given in Table I. We have calculated exponent h*(q) for other Inventions and Sinfonias 
as well, and all of them are in the 1.7 — 1.8 range (Fig. [5]). Figure . H(c) shows the width of 
singularity spectrum, /*(«*), i.e. Aa = a*(q max ) — a*(q min ) for the series is approximately, 



1.70. It value shows the power of multifractality of interevents is very strong 



41] . The values 
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of derived quantities from MF-DFAl method, are given in Table I. There are some specific 
lengths in music that have important effects in music's statistical parameters, a) Rhythm, 
b) Scale is a set of tones each having a definite pitch (perceived fundamental frequency of a 
sound) and each having a specific frequency ratio compared to the others, that is, each having 
specific interval relative to all the other pitches, arranged in a sequence from low to high, 
or alternatively, from high to low. Composers often transform musical patterns by moving 
every note in the pattern by a constant number of scale intervals. Since the intervals of a 
scale can have various sizes, this process introduces subtle melodic and harmonic variation 
into the music. This variation is what gives scalar music much of its complexity), c) Bar or 
measure is a segment of time defined as a given number of beats (the basic time unit of a 
piece) of a given duration. More over in the Persian music (Eastern music), there are others 
characteristic lengths such as Gushe which is based on the change of the tonic or stop notes. 
Existence of these characteristic lengths can effect on music complexity. The music style 
and more importantly composers can use various the characteristic lengths in the music. We 
check a few pieces of other genders of music like jazz and Persian traditional music (Eastern 
music) that obtained h*(2) in ranges, (1.75 — 1.9) and (1.3 — 1.7), respectively. We find more 
characteristic lengths in the Persian traditional music which leads to decrease the exponents. 

Usually, in MF-DFA method, deviation from a straight line in log- log plot of eq. ([6]) occurs 
for small scales s. This deviation limits capability of DFA to determine correct correlation 
behavior for very short scales and in regime of small s. The modified MF-DFA is defined as 
follows Ul: 



^ _ jF shuf (s0]2)1 /2 s l/2 



(W»l), 



(16) 



where ([i^ huf (s)] 2 ) 1//2 denotes the usual MF-DFA fluctuation function, defined in eq.flHD, 
averaged over several configurations of shuffled data taken from original time series, and 
s' ~ iV/40. The values of h*(q) exponent obtained by modified MF-DFAl methods for 
frequency series time series is 1.70 ±0.03. The relative deviation of the new exponent which 
is obtained by modified MF-DFAl in comparison to MF-DFAl for original data is less than 
5%. 

Now, we are interested in to determine source of multifractality. In general, two different 
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FIG. 6: The new exponent, h*(q), after double profiling, as a function of q for the original, 
surrogate and shuffled data. 

types of multifractality in time series can be distinguished: (i) multifractality due to a 
fatness of probability density function (PDF) of the time series. In this case multifractality 
cannot be removed by shuffling series, (ii) multifractality due to different correlations in 
small and large scale fluctuations. In this case data may have a PDF with finite moments, 
e.g. a Gaussian distribution. Thus corresponding shuffled time series will exhibit mono- 
fractal scaling, since all long-range correlations are destroyed by shuffling procedure. If 
both kinds of multifractality are present, shuffled series will show weaker multifractality 
than original series. The easiest way to clarify the type of multifractality, is by analyzing 
corresponding shuffled and surrogate time series. Shuffling of time series destroys long-range 
correlation, therefore if multifractality belongs only to long-range correlation, we should find 
Khuiil) = 1-5. The multifractality nature due to fatness of PDF signals is not affected by 
shuffling procedure. On the other hand, to determine multifractality due to broadness of 
PDF, phase of discrete fourier transform (DFT) coefficients of frequency series time series are 
replaced with a set of pseudo independent distributed uniform (— vr, it) quantities in surrogate 
method. The correlations in surrogate series do not change, but probability function changes 
to Gaussian distribution. If multifractality in time series is due to a broad PDF, h*(q) 
obtained by surrogate method will be independent of q. If both kinds of multifractality 
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are present in frequency series time series, shuffled and surrogate series will show weaker 
multifractality than original one. 

To check nature of multifractality, we compare fluctuation function F q (s), for origi- 
nal series (after cancelation of sinusoidal trend) with result from corresponding shuffled, 
Fq hui (s) and surrogate series F| ur (s). Differences between these two fluctuation functions 
with original one, directly indicate presence of long-range correlations or broadness of prob- 
ability density function in original series. These differences can be observed in ratio plot 
of F q (s) I Fq hnf (s) and F q (s)/F* m (s) with respect s [40j]. Since anomalous scaling due to a 
broad probability density affects both of F q (s) and F q hui (s) in same way, only multifractality 
due to correlations will be observed in F q (s) / F* hui (s) . Scaling behavior of these ratios are 

F ? (s)/Ff uf (s) ~ s^W-^hufW = S ^ OT (?) ) (17) 

F q (s)/F» m (s) ~ s h *(i)- h suM = s h* PDF (i). (18) 

If only fatness of PDF is responsible for multifractality, one should have h*(q) = hl huf (q) 
and hl OT (q) = 0. On the other hand, deviations from h* or (q) = indicates existence of 
correlations, and q dependence of h* or (q) indicates that multifractality is due to long rage 
correlation. If only correlation multifractality is present, one finds h* hu{ (q) = 1-5- If both 
distribution and correlation multifractality are present, both of h* hui (q) and hl ur (q) will 
depend on q. The q dependence of exponent h*(q) for original, surrogate and shuffled time 
series are shown in Figures [61 q dependence of h COT and h PDF shows that multifractality 
nature of frequency series time series is due to both broadness of PDF and long-range 
correlation. Absolute value of h PDF (q) is greater than h* OT (q), so multifractality due to 
correlation is weaker than mulifractality due to fatness. 

Deviation of hl UT (q) and /t* huf (q) from h*(q) can be determined by using x 2 test as follows: 



N 



symbol " o " can be replaced by "sur" and "shuf , to determine the confidence level of ^i* ur 
and /i* huf to the new exponents, h*(q), of original series, respectively. Reduced chi-square 

2 

X 2 vo = (A/" is number of degree of freedom) for shuffled and surrogate time series are 
16.07, 21.30, respectively. Width of singularity spectrum, /*(«*), i.e. Act = a*(q m i n ) — 
ot* (imax) for original, surrogate and shuffled time series are approximately, 1.04, 0.65 and 



0.92 respectively. These values conclude that multifractality due to fat tail is dominant 
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TABLE I: Values of h*(q = 2), r*(q = 2) exponents and Width of singularity spectrum, f*(a*), 
Aa for q = 2.0 of Invention no.l obtained by MF-DFAl. 





h*(2) 


T*(2) 


Aa 


Original 


1.74 ±0.03 


2.48 ± 0.03 


1.70 


Surrogate 


1.71 ±0.03 


2.44 ±0.03 


0.80 


Shuffled 


0.50 ±0.03 


0.02 ±0.03 


0.50 



Values of the exponents, h*(q = 2.0), r*(q = 2) and width of singularity spectrum, 
f*(a*), Aa for the original, shuffled and surrogate of frequency series obtained with MF- 
DFAl method are reported in Table I. 

IV. CONCLUSION 

MF-DFA method allows us to determine multifractal characterization of non- stationary 
and stationary time series. We have shown that MF-DFAl result of frequency series of 
Bach Inventions. Applying MF-DFAl method is demonstrated that frequency series have 
long term correlation. We calculated the second moment exponent after the new profile for 
other Inventions and Sinfonias and they are in the 1.7 — 1.8 range, q dependence of h*(q) 
and T*(q), shows that frequency series has multifractal behavior. By comparing the second 
moment exponent of original time series with shuffled and surrogate one's, we have found 
that multifractality due to broadness of probability density function has more contribution 
than correlation in Inventions. 
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